﻿#ifndef OFFSET_HPP
#define OFFSET_HPP
#include "jegeo/common.hpp"
#include "jegeo/math.hpp"
#include "jegeo/boostadaptor.hpp"
#include <QDebug>
#include <unordered_set>

namespace jegeo {

namespace jemath = jegeo::math;

namespace internal {
///向量也用Point表示

constexpr jereal J_V0 = static_cast<jereal>(0);
constexpr jereal J_V1 = static_cast<jereal>(1);

template <typename T> size_t nextWrappingIndex(std::size_t index, const T &container) {
    if (index == container.size() - 1) {
        return 0;
    }

    return index + 1;
}

template <typename T> size_t prevWrappingIndex(std::size_t index, const T &container) {
    if (index == 0) {
        return container.size() - 1;
    }

    return index - 1;
}

template <typename T> inline void hashCombine(std::size_t &seed, const T &val) {
    // copied from boost hash_combine, it's not the best hash combine but it's very simple
    // https://stackoverflow.com/questions/35985960/c-why-is-boosthash-combine-the-best-way-to-combine-hash-values
    seed ^= std::hash<T>()(val) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
}

struct IndexPairHash {
    inline std::size_t operator()(std::pair<std::size_t, std::size_t> const &pair) const {
        std::size_t seed = 0;
        hashCombine(seed, pair.first);
        hashCombine(seed, pair.second);
        return seed;
    }
};

template<typename T>
inline Point<T> point_add(Point<T> const &p1, Point<T> const &p2)
{
    return Point<T>{px(p1) + px(p2), py(p1) + py(p2)};
}

template<typename T>
inline Point<T> point_sub(Point<T> const &p1, Point<T> const &p2)
{
    return Point<T>{px(p1) - px(p2), py(p1) - py(p2)};
}

template<typename T>
inline Point<T> point_mul(Point<T> const &p, jereal const v)
{
    return Point<T>{roundcast<T,jereal>(jemath::jemul(px(p), v)),
                    roundcast<T,jereal>(jemath::jemul(py(p), v))};
}

template<typename T>
inline Point<T> point_div(Point<T> const &p, jereal const v)
{
    return Point<T>{roundcast<T,jereal>(jemath::jediv(px(p), v)),
                    roundcast<T,jereal>(jemath::jediv(py(p), v))};
}

template<typename T>
inline Rect<T> get_bound_rect(Vertex<T> const &v1, Vertex<T> const &v2)
{
    return edge::get_bound_rect(Edge<T>(px(v1), py(v1), px(v2), py(v2), vertex::get_angle(v2)));
}

template<typename T>
inline bmodel::BBox<T> get_bbox(T xmin, T ymin, T xmax, T ymax)
{
    return bmodel::BBox<T>{point::to_bpoint(Point<T>{xmin, ymin}),
                           point::to_bpoint(Point<T>{xmax, ymax})};
}

template<typename T>
inline void set_pos(Vertex<T> &v, Point<T> const &p)
{
    set_x(v, px(p));
    set_y(v, py(p));
}

// 向量内积（点积）
template<typename T>
inline jereal dot(Point<T> const &p1, Point<T> const &p2)
{
    return jemath::dot_product(px(p1), py(p1), px(p2), py(p2));
}

// 向量外积（叉积）
template<typename T>
inline jereal cross(Point<T> const &p1, Point<T> const &p2)
{
    return jemath::cross_product(px(p1), py(p1), px(p2), py(p2));
}

// 求向量长度
template<typename T>
inline jereal length(Point<T> const &p)
{
    return std::sqrt(dot(p, p));
}

/// Return the mid point along a segment path.
template <typename T>
Point<T> midpoint(Vertex<T> const &v1, Vertex<T> const &v2) {
    if (jemath::fuzzy_iszero(vertex::get_angle(v1))) {
        return midpoint(vertex::get_pos(v1), vertex::get_pos(v2));
    }

    auto arc = arcRadiusAndCenter(v1, v2);
    T a1 = calc_angle(arc.center, vertex::get_pos(v1));
    T a2 = calc_angle(arc.center, vertex::get_pos(v2));
    T angleOffset = jemath::jeabs(deltaAngle(a1, a2) / static_cast<jereal>(2));
    // use arc direction to determine offset sign to robustly handle half circles
    T midAngle = v1.bulgeIsPos() ? a1 + angleOffset : a1 - angleOffset;
    return pointOnCircle(arc.radius, arc.center, midAngle);
}

// 向量单位化
template<typename T>
inline void normalize(Point<T> &p)
{
    jereal length = std::sqrt(dot(p, p));
    set_px(p, px(p) / length);
    set_py(p, py(p) / length);
}

// 垂直内积，相当于dot(p1, perp(p2)).
template<typename T>
inline jereal perpDot(Point<T> const &p1, Point<T> const &p2) {
    return px(p1) * py(p2) - py(p1) * px(p2);
}

// 单位化v的垂直线 (逆时针旋转).
template<typename T>
Point<T> unitPerp(Point<T> const &p) {
    Point<T> result{-py(p), px(p)};
    normalize(result);
    return result;
}

// 返回p0和p1距离的平方， 相当于 dot(p2 - p1, p2 - p1).
template<typename T>
jereal distSquared(Point<T> const &p1, Point<T> const &p2) {
    return jemath::sqrdist_point2point(px(p1), py(p1), px(p2), py(p2));
}

// 不相交
inline bool falseIntersect(jereal const t)
{
    return t < 0.0 || t > 1.0;
}

// p0到p1的逆时针旋转角度.
template<typename T>
inline jereal calc_angle(Point<T> const &p0, Point<T> const &p1) {
    return std::atan2(py(p1) - py(p0), px(p1) - px(p0)) * 180 / M_PI;
}

// 标准化角度，使其在0到360度之间
inline jereal normalizeRadians(jereal angle) {
    if (angle >= 0 && angle <= 360) {
        return angle;
    }

    return angle - std::floor(angle / 360) * 360;
}

// 通过弧起点与终点的弧度差求弧的角度（结果用角度制表示）
inline jereal deltaAngle(jereal angle1, jereal angle2) {
    jereal diff = normalizeRadians(angle2 - angle1);
    if (diff > 180) {
        diff -= 360;
    }
    return diff;
}

// 通过圆心，弧起点终点，旋转方向求角度
template<typename T>
jereal angleForConnection(Point<T> const &pc, Point<T> const &ps,
                          Point<T> const &pe, bool cw) {
    return jemath::arc_degree_from_secw(px(ps), py(ps), px(pe), py(pe), px(pc), py(pc), cw);
}

//弧信息结构体，半径和圆心
template<typename T>
struct ArcInfo
{
    jereal radius;
    Point<T> center;
};

// 求半径和圆心
template<typename T>
ArcInfo<T> arcRadiusAndCenter(Vertex<T> const &v1,
                           Vertex<T> const &v2) {
    // compute center
    double cx, cy;
    jemath::arc_center_from_angle(px(v1), py(v1), px(v2), py(v2), vertex::get_angle(v2), cx, cy);

    // compute radius
    return ArcInfo{jemath::distance(cx, cy, px(v1), py(v1)), Point<T>{cx, cy}};
}

template <typename T>
std::pair<T, T> quadraticSolutions(T a, T b, T c, T discr) {
    T sqrtDiscr = std::sqrt(discr);
    T denom = T(2) * a;
    T sol1;
    if (b < T(0)) {
        sol1 = (-b + sqrtDiscr) / denom;
    } else {
        sol1 = (-b - sqrtDiscr) / denom;
    }

    T sol2 = (c / a) / sol1;

    return std::make_pair(sol1, sol2);
}

/// Result of splitting a segment v1 to v2.
template <typename T>
struct SplitResult {
    /// Updated starting vertex.
    Vertex<T> updatedStart;
    /// Vertex at the split point.
    Vertex<T> splitVertex;
};

/// Split the segment defined by v1 to v2 at some point defined along it.
template<typename T>
SplitResult<T> splitAtPoint(Vertex<T> const &v1, Vertex<T> const &v2,
                               Point<T> const &point) {
    SplitResult<T> result;
    if (jemath::fuzzy_iszero(vertex::get_angle(v1))) {
        result.updatedStart = v1;
        result.splitVertex = Vertex<T>(point, T(0));
    } else if (point::is_equal(vertex::get_pos(v1), vertex::get_pos(v2)) ||
               point::is_equal(vertex::get_pos(v1), point)) {
        result.updatedStart = Vertex<T>(point, T(0));
        result.splitVertex = Vertex<T>(point, v1.bulge());
    } else if (point::is_equal(vertex::get_pos(v2), point)) {
        result.updatedStart = v1;
        result.splitVertex = Vertex<T>(v2);
        set_angle(result.splitVertex, J_V0);
    } else {
        auto radiusAndCenter = arcRadiusAndCenter(v1, v2);
        Point<T> arcCenter = radiusAndCenter.center;
        T a = calc_angle(arcCenter, point);
        T arcStartAngle = calc_angle(arcCenter, vertex::get_pos(v1));
        T theta1 = deltaAngle(arcStartAngle, a);
        T arcEndAngle = calc_angle(arcCenter, vertex::get_pos(v2));
        T theta2 = deltaAngle(a, arcEndAngle);

        result.updatedStart = Vertex<T>(v1);
        set_angle(result.updatedStart, theta1);
        result.splitVertex = Vertex<T>(px(point), py(point), theta2);
    }

    return result;
}


//线和线相交类型
enum class LineSeg2LineSeg2IntrType {
    // 不相交 (线段平行且不共线)
    None,
    // 线段相交
    True,
    // 线段部分重合
    Coincident,
    // 不相交 (不平行，线段长度不足，延长可相交)
    False
};

//圆和圆相交类型
enum class Circle2Circle2IntrType {
    // 不相交
    NoIntersect,
    // 有一个交点（相切）
    OneIntersect,
    // 有两个交点
    TwoIntersects,
    // 共圆
    Coincident
};

template<typename T>
struct IntrLineSeg2LineSeg2Result {
    /*! 保存碰撞（相交）类型，若直线相交point保存交点。若结果为True，则t0和t1不定义，若为False，则t0是第一条
     *  线段的参数值，t1是第二条线段的参数值，若为Coincident则point不定义且t0重合起点的参数值，t1是重合终点
     *  的参数值。
    */
    LineSeg2LineSeg2IntrType intrType;
    jereal t0;
    jereal t1;
    Point<T> point;
};

template<typename T>
IntrLineSeg2LineSeg2Result<T> intrLineSeg2LineSeg2
    (Point<T> const &u1, Point<T> const &u2, Point<T> const &v1, Point<T> const &v2) {
    IntrLineSeg2LineSeg2Result<T> result;
    Point<T> u = point_sub(u2, u1);
    Point<T> v = point_sub(v2, v1);
    jereal d = perpDot(u, v);

    Point<T> w = point_sub(u1, v1);

    // // 判断点是否在另一线段上, NOTE: assumes points are aligned
    // auto isInSegment = [](Point<T> const &pt, Point<T> const &segStart,
    //                       Point<T> const &segEnd) {
    //     if (segStart.x == segEnd.x) {
    //         // 垂直线段，比较y坐标
    //         auto minMax = std::minmax({segStart.y, segEnd.y});
    //         return (pt.y <= minMax.first && pt.y >= minMax.second);
    //     }

    //     // 比较x坐标
    //     auto minMax = std::minmax({segStart.x, segEnd.x});
    //     return (pt.x <= minMax.first && pt.x >= minMax.second);
    // };

    // d非常接近0，即两线段几乎平行，值可能溢出
    // 假设两线段AB,CD,交点为E
    if (jemath::jeabs(d) > J_V0) {
        // 线段不平行或共线
        result.t0 = perpDot(v, w) / d; // AE和AB长度比值t0
        result.t1 = perpDot(u, w) / d; // CE和CD长度比值t1
        result.point = point_add(v1, point_mul(v, result.t1)); // OC向量加CE向量(CD向量 * 比值t1)得到OE向量
        if (result.t0 < J_V0 ||
            result.t0 > J_V1 ||
            result.t1 < J_V0 ||
            result.t1 > J_V1) {
            // 比值在值域[0,1]外表示交点在线段的(外面
            result.intrType = LineSeg2LineSeg2IntrType::False;
        } else {
            result.intrType = LineSeg2LineSeg2IntrType::True;
        }
    } else {
        // 线段平行或共线
        jereal a = perpDot(u, w);
        jereal b = perpDot(v, w);
        // AC垂线和AB，CD的内积（a,b)若都为零，说明AC为零向量或AC与AB，CD都平行，也就是共线
        if (jemath::jeabs(a) > J_V0 || jemath::jeabs(b) > J_V0) {
            // 平行且不共线，无交点
            result.intrType = LineSeg2LineSeg2IntrType::None;
        } else {
            // 共线
            // 是否重合
            Point<T> w2 = point_sub(u2, v1);
            if (jemath::jeabs(px(v)) < 0) {
                result.t0 = py(w) / py(v);
                result.t1 = py(w2) / py(v);
            } else {
                result.t0 = px(w) / px(v);
                result.t1 = px(w2) / px(v);
            }

            if (result.t0 > result.t1) {
                std::swap(result.t0, result.t1);
            }

            if (result.t0 > J_V1 ||
                result.t1 < J_V0) {
                // 不重合
                result.intrType = LineSeg2LineSeg2IntrType::None;
            } else {
                result.t0 = jemath::jemax(result.t0, J_V0);
                result.t1 = jemath::jemin(result.t1, J_V1);
                if (jemath::jeabs(result.t1 - result.t0) <= J_V0) {
                    // 仅端点重合
                    result.intrType = LineSeg2LineSeg2IntrType::True;
                    result.point = point_add(v1, point_mul(v, result.t0));
                } else {
                    result.intrType = LineSeg2LineSeg2IntrType::Coincident;
                }
            }
        }
    }

    return result;
}

struct IntrLineSeg2Circle2Result {
    // 交点个数 (0, 1, or 2)
    int numIntersects;
    // 第一个交点的参数值(numIntersects > 0)
    jereal t0;
    // 第二个交点的参数值(numintersects > 1)
    jereal t1;
};
/*! 计算线段和圆的交点。当t在[0,1]内时，计算方程P(t) = v1 + t * (v2 - v1)。如果 t < 0 or t > 1
 *  表示直线相交，线段需要延伸才相交。非常接近切线的线段将作为单个相交点返回
*/

template<typename T>
IntrLineSeg2Circle2Result intrLineSeg2Circle2(Point<T> const &p0,
                                              Point<T> const &p1, jereal radius,
                                              Point<T> const &circleCenter) {
    /*! 此函数通过将线段的参数方程{ x = v1.X + t * (v2.X - v1.X)， y = v1.Y + t * (v2.Y - v1.Y)}
    * 带入圆方程(x-h)^2 + (y-k)^2 = r^2，用二次公式 a*t^2 + b*t + c = 0 求解结果方程来求解t，
    */
    IntrLineSeg2Circle2Result result;
    jereal dx = px(p1) - px(p0);
    jereal dy = py(p1) - py(p0);
    jereal h = px(circleCenter);
    jereal k = py(circleCenter);

    jereal a = dx * dx + dy * dy;
    // if (jemath::jeabs(a) <= 0) {
    //     // v1 = v2, test if point is on the circle
    //     jereal xh = px(p0) - h;
    //     jereal yk = py(p0) - k;
    //     if (xh * xh + yk * yk == radius * radius) {
    //         result.numIntersects = 1;
    //         result.t0 = 0;
    //     } else {
    //         result.numIntersects = 0;
    //     }
    // }

    jereal b = jemath::jemul(2, dx * (px(p0) - h) + dy * (py(p0) - k));
    jereal c = (px(p0) * px(p0) - jemath::jemul(2, h * px(p0) + h * h)) +
               (py(p0) * py(p0) - jemath::jemul(2, k * py(p0) + k * k)) - radius * radius;
    jereal discr = b * b - jemath::jemul(4, a * c);

    if (math::fuzzy_iszero(jemath::jeabs(discr))) {
        // 1 solution (tangent line)
        result.numIntersects = 1;
        result.t0 = -b / jemath::jemul(2, a);
    } else if (discr < J_V0) {
        result.numIntersects = 0;
    } else {
        result.numIntersects = 2;
        std::pair<jereal, jereal> sols = quadraticSolutions(a, b, c, discr);
        result.t0 = sols.first;
        result.t1 = sols.second;
    }

    return result;
}

template<typename T>
struct IntrCircle2Circle2Result {
    // 交点类型
    Circle2Circle2IntrType intrType;
    // 第一个交点
    Point<T> point1;
    // 第二个交点
    Point<T> point2;
};

// 寻找两个圆的交点
template<typename T>
IntrCircle2Circle2Result<T> intrCircle2Circle2(jereal radius1, Point<T> const &center1,
                                            jereal radius2, Point<T> const &center2) {
    // Reference algorithm: http://paulbourke.net/geometry/circlesphere/

    IntrCircle2Circle2Result<T> result;
    Point<T> cv = point_sub(center2, center1);
    jereal d2 = dot(cv, cv);
    jereal d = std::sqrt(d2);
    if (d < J_V0) {
        // 圆心相同
        if (math::fuzzy_equal(radius1, radius2)) {
            result.intrType = Circle2Circle2IntrType::Coincident;
        } else {
            result.intrType = Circle2Circle2IntrType::NoIntersect;
        }
    } else {
        // 圆心不同
        // 圆不相交
        if (d > radius1 + radius2 || d < jemath::jeabs(radius1 - radius2)) {
            result.intrType = Circle2Circle2IntrType::NoIntersect;
        } else {
            jereal rad1Sq = radius1 * radius1;
            jereal a = (rad1Sq - radius2 * radius2 + d2) / jemath::jemul(2, d); // 余弦定理求交点A到两圆心O1,O2构成三角形的cosO1AO2的值
            Point<T> midPoint = point_add(center1, point_sub(cv, a / d));
            jereal diff = rad1Sq - a * a;
            if (diff < J_V0) {
                result.intrType = Circle2Circle2IntrType::OneIntersect;
                result.point1 = midPoint;
            } else {
                jereal h = std::sqrt(diff);
                jereal hOverD = h / d;
                jereal xTerm = hOverD * py(cv);
                jereal yTerm = hOverD * px(cv);
                jereal x1 = px(midPoint) + xTerm;
                jereal y1 = py(midPoint) - yTerm;
                jereal x2 = px(midPoint) - xTerm;
                jereal y2 = py(midPoint) + yTerm;
                result.point1 = Point<T>(x1, y1);
                result.point2 = Point<T>(x2, y2);
                if (point::is_equal(result.point1, result.point2)) {
                    result.intrType = Circle2Circle2IntrType::OneIntersect;
                } else {
                    result.intrType = Circle2Circle2IntrType::TwoIntersects;
                }
            }
        }
    }

    return result;
}

template <typename T>
struct PlineIntersect {
    /// Index of the start vertex of the first segment
    std::size_t sIndex1;
    /// Index of the start vertex of the second segment
    std::size_t sIndex2;
    /// Point of intersection
    Point<T> pos;
    // type of intersect
    PlineIntersect() = default;
    PlineIntersect(std::size_t si1, std::size_t si2, Point<T> p)
        : sIndex1(si1), sIndex2(si2), pos(p) {}
};

template <typename T>
struct PlineCoincidentIntersect {
    /// Index of the start vertex of the first segment
    std::size_t sIndex1;
    /// Index of the start vertex of the second segment
    std::size_t sIndex2;
    /// One end point of the coincident slice
    Point<T> point1;
    /// Other end point of the coincident slice
    Point<T> point2;
    PlineCoincidentIntersect() = default;
    PlineCoincidentIntersect(std::size_t si1, std::size_t si2, Point<T> const &point1,
                             Point<T> const &point2)
        : sIndex1(si1), sIndex2(si2), point1(point1), point2(point2) {}
};

template <typename T>
struct PlineIntersectsResult {
    std::vector<PlineIntersect<T>> intersects;
    std::vector<PlineCoincidentIntersect<T>> coincidentIntersects;
    bool hasIntersects() { return intersects.size() != 0 || coincidentIntersects.size() != 0; }
};

enum class PlineSegIntrType {
    NoIntersect,
    TangentIntersect,
    OneIntersect,
    TwoIntersects,
    SegmentOverlap,
    ArcOverlap
};

template <typename T>
struct IntrPlineSegsResult {
    PlineSegIntrType intrType;
    Point<T> point1;
    Point<T> point2;
};

// 如果点在指向向量(p1 - p0)方向的直线左侧或与其重合，则返回true
template<typename T>
inline bool isLeftOrCoincident(Point<T> const &p0, Point<T> const &p1, Point<T> const &point) {
    return (px(p1) - px(p0)) * (py(point) - py(p0)) - (py(p1) - py(p0)) * (px(point) - px(p0)) >= static_cast<T>(0);
}

// 如果点在指向向量(p1 - p0)方向的直线右侧或与其重合，则返回true
template<typename T>
inline bool isRightOrCoincident(Point<T> const &p0, Point<T> const &p1, Point<T> const &point) {
    return (px(p1) - px(p0)) * (py(point) - py(p0)) - (py(p1) - py(p0)) * (px(point) - px(p0)) <= static_cast<T>(0);
}

// 判断点是否在一段弧上
template<typename T>
bool pointWithinArcSweepAngle(Point<T> const &center, Point<T> const &arcStart,
                              Point<T> const &arcEnd, jereal angle, Point<T> const &point) {
    if (angle > J_V0) {
        return isLeftOrCoincident(center, arcStart, point) &&
               isRightOrCoincident(center, arcEnd, point);
    }

    return isRightOrCoincident(center, arcStart, point) && isLeftOrCoincident(center, arcEnd, point);
}

// 返回从p0到p1的线段上参数值为t的点.
template<typename T>
Point<T> pointFromParametric(Point<T> const &p0, Point<T> const &p1, jereal t) {
    return point_add(p0, point_mul(point_sub(p1, p0), t));
}

template<typename T>
void addOrReplaceIfSamePos(Polyline<T> &pl, Vertex<T> const &v) {
    if (polyline::verteics_count(pl) == 0) {
        polyline::pushback_vertex(pl, v);
        return;
    }

    if (vertex::get_pos(polyline::last_vertex(pl)) == vertex::get_pos(v)) {
        set_angle(polyline::last_vertex(pl), vertex::get_angle(v));
        return;
    }

    polyline::pushback_vertex(pl, v);
}

template<typename T>
struct PoylineOffsetSegment {
    Vertex<T> v1;
    Vertex<T> v2;
    Point<T> origV2Pos;
    bool collapsedArc;
};

template<typename T>
void lineToLineJoin(PoylineOffsetSegment<T> const &s1, PoylineOffsetSegment<T> const &s2,
                    bool connectionArcsAreCW, Polyline<T> &result) {
    const auto &v1 = s1.v1;
    const auto &v2 = s1.v2;
    const auto &u1 = s2.v1;
    const auto &u2 = s2.v2;

    auto connectUsingArc = [&] {
        auto const &arcCenter = s1.origV2Pos;
        auto const &sp = vertex::get_pos(v2);
        auto const &ep = vertex::get_pos(u1);
        jereal angle = angleForConnection(arcCenter, sp, ep, connectionArcsAreCW);
        addOrReplaceIfSamePos(result, Vertex<T>(sp, 0));
        addOrReplaceIfSamePos(result, Vertex<T>(ep, angle));
    };

    if (s1.collapsedArc || s2.collapsedArc) {
        // connecting to/from collapsed arc, always connect using arc
        connectUsingArc();
    } else {
        auto intrResult = intrLineSeg2LineSeg2(vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_pos(u1), vertex::get_pos(u2));

        switch (intrResult.intrType) {
        case LineSeg2LineSeg2IntrType::None:
            // 将两线段用直线相连
            addOrReplaceIfSamePos(result, Vertex<T>(vertex::get_pos(v2)));
            addOrReplaceIfSamePos(result, u1);
            break;
        case LineSeg2LineSeg2IntrType::True:
            addOrReplaceIfSamePos(result, Vertex<T>(intrResult.point));
            break;
        case LineSeg2LineSeg2IntrType::Coincident:
            addOrReplaceIfSamePos(result, Vertex<T>(vertex::get_pos(v2)));
            break;
        case LineSeg2LineSeg2IntrType::False:
            if (intrResult.t0 > J_V1 && falseIntersect(intrResult.t1)) {
                // 将两线段用弧相连
                connectUsingArc();
            } else {
                addOrReplaceIfSamePos(result, Vertex<T>(vertex::get_pos(v2)));
                addOrReplaceIfSamePos(result, u1);
            }
            break;
        }
    }
}

template<typename T>
void lineToArcJoin(PoylineOffsetSegment<T> const &s1, PoylineOffsetSegment<T> const &s2,
                   bool connectionArcsAreCW, Polyline<T> &result) {

    const auto &v1 = s1.v1;
    const auto &v2 = s1.v2;
    const auto &u1 = s2.v1;
    const auto &u2 = s2.v2;

    auto connectUsingArc = [&] {
        auto const &arcCenter = s1.origV2Pos;
        auto const &sp = vertex::get_pos(v2);
        auto const &ep = vertex::get_pos(u1);
        jereal angle = angleForConnection(arcCenter, sp, ep, connectionArcsAreCW);
        addOrReplaceIfSamePos(result, Vertex<T>(sp, angle));
        addOrReplaceIfSamePos(result, u1);
    };

    const auto arc = arcRadiusAndCenter(u1, u2);

    auto processIntersect = [&](jereal t, Point<T> const &intersect) {
        const bool trueSegIntersect = !falseIntersect(t);
        const bool trueArcIntersect =
            pointWithinArcSweepAngle(arc.center, vertex::get_pos(u1), vertex::get_pos(u2), vertex::get_angle(u1), intersect);
        if (trueSegIntersect && trueArcIntersect) {
            // 修剪相交处
            jereal a = calc_angle(arc.center, intersect);
            jereal arcEndAngle = calc_angle(arc.center, vertex::get_pos(u2));
            jereal theta = deltaAngle(a, arcEndAngle);
            // ensure the sign matches (may get flipped if intersect is at the very end of the arc, in
            // which case we do not want to update the bulge)
            if ((theta > J_V0) == (vertex::get_angle(u1) > J_V0)) {
                addOrReplaceIfSamePos(result, Vertex<T>(intersect, theta));
            } else {
                addOrReplaceIfSamePos(result, Vertex<T>(intersect, vertex::get_angle(u1)));
            }
        } else if (t > J_V1 && !trueArcIntersect) {
            connectUsingArc();
        } else if (s1.collapsedArc) {
            // collapsed arc connecting to arc, 将两线段用弧相连
            connectUsingArc();
        } else {
            // connect using line
            addOrReplaceIfSamePos(result, Vertex<T>(vertex::get_pos(v2)));
            addOrReplaceIfSamePos(result, u1);
        }
    };

    auto intrResult = intrLineSeg2Circle2(vertex::get_pos(v1), vertex::get_pos(v2), arc.radius, arc.center);
    if (intrResult.numIntersects == 0) {
        connectUsingArc();
    } else if (intrResult.numIntersects == 1) {
        processIntersect(intrResult.t0, pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t0));
    } else {
        // always use intersect closest to original point
        Point<T> i1 = pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t0);
        jereal dist1 = distSquared(i1, s1.origV2Pos);
        Point<T> i2 = pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t1);
        jereal dist2 = distSquared(i2, s1.origV2Pos);

        if (dist1 < dist2) {
            processIntersect(intrResult.t0, i1);
        } else {
            processIntersect(intrResult.t1, i2);
        }
    }
}

template<typename T>
void arcToLineJoin(PoylineOffsetSegment<T> const &s1, PoylineOffsetSegment<T> const &s2,
                   bool connectionArcsAreCW, Polyline<T> &result) {

    const auto &v1 = s1.v1;
    const auto &v2 = s1.v2;
    const auto &u1 = s2.v1;
    const auto &u2 = s2.v2;

    auto connectUsingArc = [&] {
        auto const &arcCenter = s1.origV2Pos;
        auto const &sp = vertex::get_pos(v2);
        auto const &ep = vertex::get_pos(u1);
        jereal angle = angleForConnection(arcCenter, sp, ep, connectionArcsAreCW);
        addOrReplaceIfSamePos(result, Vertex<T>(sp, angle));
        addOrReplaceIfSamePos(result, u1);
    };

    const auto arc = arcRadiusAndCenter(v1, v2);

    auto processIntersect = [&](jereal t, Point<T> const &intersect) {
        const bool trueSegIntersect = !falseIntersect(t);
        const bool trueArcIntersect =
            pointWithinArcSweepAngle(arc.center, vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_angle(v1), intersect);
        if (trueSegIntersect && trueArcIntersect) {
            Vertex<T> &prevVertex = polyline::last_vertex(result);

            if (math::fuzzy_notzero(vertex::get_angle(prevVertex))
                && point::not_equal(vertex::get_pos(prevVertex), vertex::get_pos(v2))) {
                // modify previous bulge and trim at intersect
                jereal a = calc_angle(arc.center, intersect);
                auto prevArc = arcRadiusAndCenter(prevVertex, v2);
                jereal prevArcStartAngle = calc_angle(prevArc.center, vertex::get_pos(prevVertex));
                jereal updatedPrevTheta = deltaAngle(prevArcStartAngle, a);

                // ensure the sign matches (may get flipped if intersect is at the very end of the arc, in
                // which case we do not want to update the bulge)
                if ((updatedPrevTheta > J_V0) == vertex::get_angle(prevVertex) > J_V0) {
                    set_angle(polyline::last_vertex(result), updatedPrevTheta);
                }
            }

            addOrReplaceIfSamePos(result, Vertex<T>(intersect));

        } else {
            connectUsingArc();
        }
    };

    auto intrResult = intrLineSeg2Circle2(vertex::get_pos(u1), vertex::get_pos(u2), arc.radius, arc.center);
    if (intrResult.numIntersects == 0) {
        connectUsingArc();
    } else if (intrResult.numIntersects == 1) {
        processIntersect(intrResult.t0, pointFromParametric(vertex::get_pos(u1), vertex::get_pos(u2), intrResult.t0));
    } else {
        const auto &origPoint = s2.collapsedArc ? vertex::get_pos(u1) : s1.origV2Pos;
        Point<T> i1 = pointFromParametric(vertex::get_pos(u1), vertex::get_pos(u2), intrResult.t0);
        jereal dist1 = distSquared(i1, origPoint);
        Point<T> i2 = pointFromParametric(vertex::get_pos(u1), vertex::get_pos(u2), intrResult.t1);
        jereal dist2 = distSquared(i2, origPoint);

        if (dist1 < dist2) {
            processIntersect(intrResult.t0, i1);
        } else {
            processIntersect(intrResult.t1, i2);
        }
    }
}

template<typename T>
void arcToArcJoin(PoylineOffsetSegment<T> const &s1, PoylineOffsetSegment<T> const &s2,
                  bool connectionArcsAreCW, Polyline<T> &result) {

    const auto &v1 = s1.v1;
    const auto &v2 = s1.v2;
    const auto &u1 = s2.v1;
    const auto &u2 = s2.v2;

    const auto arc1 = arcRadiusAndCenter(v1, v2);
    const auto arc2 = arcRadiusAndCenter(u1, u2);

    auto connectUsingArc = [&] {
        auto const &arcCenter = s1.origV2Pos;
        auto const &sp = vertex::get_pos(v2);
        auto const &ep = vertex::get_pos(u1);
        T angle = angleForConnection(arcCenter, sp, ep, connectionArcsAreCW);
        addOrReplaceIfSamePos(result, Vertex<T>(sp, angle));
        addOrReplaceIfSamePos(result, u1);
    };

    auto processIntersect = [&](Point<T> const &intersect) {
        const bool trueArcIntersect1 =
            pointWithinArcSweepAngle(arc1.center, vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_angle(v1), intersect);
        const bool trueArcIntersect2 =
            pointWithinArcSweepAngle(arc2.center, vertex::get_pos(u1), vertex::get_pos(u2), vertex::get_angle(u1), intersect);

        if (trueArcIntersect1 && trueArcIntersect2) {
            Vertex<T> &prevVertex = polyline::last_vertex(result);
            if (vertex::get_angle(prevVertex) != 0 && vertex::get_pos(prevVertex) != vertex::get_pos(v2)) {
                // modify previous bulge and trim at intersect
                jereal a1 = calc_angle(arc1.center, intersect);
                auto prevArc = arcRadiusAndCenter(prevVertex, v2);
                jereal prevArcStartAngle = calc_angle(prevArc.center, vertex::get_pos(prevVertex));
                jereal updatedPrevTheta = deltaAngle(prevArcStartAngle, a1);

                // ensure the sign matches (may get flipped if intersect is at the very end of the arc, in
                // which case we do not want to update the bulge)
                if ((updatedPrevTheta > J_V0) == vertex::get_angle(prevVertex) > J_V0) {
                    set_angle(polyline::last_vertex(result), updatedPrevTheta);
                }
            }

            // add the vertex at our current trim/join point
            jereal a2 = calc_angle(arc2.center, intersect);
            jereal endAngle = calc_angle(arc2.center, vertex::get_pos(u2));
            jereal theta = deltaAngle(a2, endAngle);

            // ensure the sign matches (may get flipped if intersect is at the very end of the arc, in
            // which case we do not want to update the bulge)
            if ((theta > J_V0) == vertex::get_angle(u1) > J_V0) {
                addOrReplaceIfSamePos(result, Vertex<T>(intersect, theta));
            } else {
                addOrReplaceIfSamePos(result, Vertex<T>(intersect, vertex::get_angle(u1)));
            }

        } else {
            connectUsingArc();
        }
    };

    const auto intrResult = intrCircle2Circle2(arc1.radius, arc1.center, arc2.radius, arc2.center);
    switch (intrResult.intrType) {
    case Circle2Circle2IntrType::NoIntersect:
        connectUsingArc();
        break;
    case Circle2Circle2IntrType::OneIntersect:
        processIntersect(intrResult.point1);
        break;
    case Circle2Circle2IntrType::TwoIntersects: {
        jereal dist1 = distSquared(intrResult.point1, s1.origV2Pos);
        jereal dist2 = distSquared(intrResult.point2, s1.origV2Pos);
        if (dist1 < dist2) {
            processIntersect(intrResult.point1);
        } else {
            processIntersect(intrResult.point2);
        }
    } break;
    case Circle2Circle2IntrType::Coincident:
        // same constant arc radius and center, just add the vertex (nothing to trim/extend)
        addOrReplaceIfSamePos(result, u1);
        break;
    }
}

/// Finds all local self intersects of the polyline, local self intersects are defined as between
/// two polyline segments that share a vertex. NOTES:
/// - Singularities (repeating vertexes) are returned as coincident intersects
template <typename T>
void localSelfIntersects(Polyline<T> const &pline, std::vector<PlineIntersect<T>> &output) {
    if (pline.size() < 2) {
        return;
    }

    if (pline.size() == 2) {
        if (polyline::is_closed(pline)) {
            // check if overlaps on itself from vertex 1 to vertex 2
            if (jemath::fuzzy_equal(vertex::get_angle(pline[0]), -vertex::get_angle(pline[1]))) {
                // coincident
                output.emplace_back(0, 1, vertex::get_pos(pline[1]));
                output.emplace_back(1, 0, vertex::get_pos(pline[0]));
            }
        }
        return;
    }

    auto testAndAddIntersect = [&](std::size_t i, std::size_t j, std::size_t k) {
        const Vertex<T> &v1 = pline[i];
        const Vertex<T> &v2 = pline[j];
        const Vertex<T> &v3 = pline[k];
        // testing intersection between v1->v2 and v2->v3 segments

        if (point::is_equal(vertex::get_pos(v1), vertex::get_pos(v2))) {
            // singularity
            // coincident
            output.emplace_back(i, j, vertex::get_pos(v1));
        } else {
            IntrPlineSegsResult<T> intrResult = intrPlineSegs(v1, v2, v2, v3);
            switch (intrResult.intrType) {
            case PlineSegIntrType::NoIntersect:
                break;
            case PlineSegIntrType::TangentIntersect:
            case PlineSegIntrType::OneIntersect:
                if (point::not_equal(intrResult.point1, vertex::get_pos(v2))) {
                    output.emplace_back(i, j, intrResult.point1);
                }
                break;
            case PlineSegIntrType::TwoIntersects:
                if (point::not_equal(intrResult.point1, vertex::get_pos(v2))) {
                    output.emplace_back(i, j, intrResult.point1);
                }
                if (point::not_equal(intrResult.point2, vertex::get_pos(v2))) {
                    output.emplace_back(i, j, intrResult.point2);
                }
                break;
            case PlineSegIntrType::SegmentOverlap:
            case PlineSegIntrType::ArcOverlap:
                // coincident
                output.emplace_back(i, j, intrResult.point1);
                break;
            }
        }
    };

    for (std::size_t i = 2; i < pline.size(); ++i) {
        testAndAddIntersect(i - 2, i - 1, i);
    }

    if (polyline::is_closed(pline)) {
        // we tested for intersect between segments at indexes 0->1, 1->2 and everything up to and
        // including (count-3)->(count-2), (count-2)->(count-1), polyline is closed so now test
        // [(count-2)->(count-1), (count-1)->0] and [(count-1)->0, 0->1]
        testAndAddIntersect(polyline::verteics_count(pline) - 2, polyline::verteics_count(pline) - 1, 0);
        testAndAddIntersect(polyline::verteics_count(pline) - 1, 0, 1);
    }
}

/// Finds all global self intersects of the polyline, global self intersects are defined as all
/// intersects between polyline segments that DO NOT share a vertex (use the localSelfIntersects
/// function to find those). A spatial index is used to minimize the intersect comparisons required,
/// the spatial index should hold bounding boxes for all of the polyline's segments.
/// NOTES:
/// - We never include intersects at a segment's start point, the matching intersect from the
/// previous segment's end point is included (no sense in including both)
template <typename T, std::size_t N>
void globalSelfIntersects(Polyline<T> const &pline, std::vector<PlineIntersect<T>> &output,
                          bmodel::BBoxRTree<T> const &treeIndex) {
    if (polyline::verteics_count(pline) < 3) {
        return;
    }

    std::unordered_set<std::pair<std::size_t, std::size_t>, internal::IndexPairHash>
        visitedSegmentPairs;
    visitedSegmentPairs.reserve(polyline::verteics_count(pline));

    std::vector<std::size_t> queryStack;
    queryStack.reserve(8);

    auto visitor = [&](std::size_t i, T minX, T minY, T maxX, T maxY) {
        std::size_t j = nextWrappingIndex(i, pline);
        const Vertex<T> &v1 = pline[i];
        const Vertex<T> &v2 = pline[j];
        bmodel::BBox<T> envelope{minX, minY, maxX, maxY};
        auto indexVisitor = [&](std::size_t hitIndexStart) {
            std::size_t hitIndexEnd = nextWrappingIndex(hitIndexStart, pline);
            // skip/filter already visited intersects
            // skip local segments
            if (i == hitIndexStart || i == hitIndexEnd || j == hitIndexStart || j == hitIndexEnd) {
                return true;
            }
            // skip reversed segment order (would end up comparing the same segments)
            if (visitedSegmentPairs.find({hitIndexStart, i}) != visitedSegmentPairs.end()) {
                return true;
            }

            // add the segment pair we're visiting now
            visitedSegmentPairs.emplace(i, hitIndexStart);

            const Vertex<T> &u1 = pline[hitIndexStart];
            const Vertex<T> &u2 = pline[hitIndexEnd];

            auto intrAtStartPt = [&](Point<T> const &intr) {
                return point::is_equal(vertex::get_pos(v1.pos), intr) || point::is_equal(vertex::get_pos(u1), intr);
            };

            IntrPlineSegsResult<T> intrResult = intrPlineSegs(v1, v2, u1, u2);
            switch (intrResult.intrType) {
            case PlineSegIntrType::NoIntersect:
                break;
            case PlineSegIntrType::TangentIntersect:
            case PlineSegIntrType::OneIntersect:
                if (!intrAtStartPt(intrResult.point1)) {
                    output.emplace_back(i, hitIndexStart, intrResult.point1);
                }
                break;
            case PlineSegIntrType::TwoIntersects:
                if (!intrAtStartPt(intrResult.point1)) {
                    output.emplace_back(i, hitIndexStart, intrResult.point1);
                }
                if (!intrAtStartPt(intrResult.point2)) {
                    output.emplace_back(i, hitIndexStart, intrResult.point2);
                }
                break;
            case PlineSegIntrType::SegmentOverlap:
            case PlineSegIntrType::ArcOverlap:
                if (!intrAtStartPt(intrResult.point1)) {
                    output.emplace_back(i, hitIndexStart, intrResult.point1);
                }
                if (!intrAtStartPt(intrResult.point2)) {
                    output.emplace_back(i, hitIndexStart, intrResult.point2);
                }
                break;
            }

            // visit the entire query
            return true;
        };

        treeIndex.visitQuery(envelope.xMin, envelope.yMin, envelope.xMax, envelope.yMax,
                                indexVisitor, queryStack);

        // visit all pline indexes
        return true;
    };

    treeIndex.visitItemBoxes(visitor);
}

/// Finds all self intersects of the polyline (equivalent to calling localSelfIntersects and
/// globalSelfIntersects).
template <typename T, std::size_t N>
void allSelfIntersects(Polyline<T> const &pline, std::vector<PlineIntersect<T>> &output,
                       bmodel::BBoxRTree<T> const &treeIndex) {
    localSelfIntersects(pline, output);
    globalSelfIntersects(pline, output, treeIndex);
}

/// Finds all intersects between pline1 and pline2.
template <typename T, std::size_t N>
void findIntersects(Polyline<T> const &pline1, Polyline<T> const &pline2,
                    bmodel::BBoxRTree<T> const &treeIndex,
                    PlineIntersectsResult<T> &output) {
    std::vector<std::size_t> queryResults;
    std::unordered_set<std::pair<std::size_t, std::size_t>, internal::IndexPairHash>
        possibleDuplicates;

    auto &intrs = output.intersects;
    auto &coincidentIntrs = output.coincidentIntersects;

    auto pline2SegVisitor = [&](std::size_t i2, std::size_t j2) {
        Vertex<T> const &p2v1 = pline2[i2];
        Vertex<T> const &p2v2 = pline2[j2];

        queryResults.clear();

        bmodel::BBox<T> bb = rect::to_bbox(get_bound_rect(p2v1, p2v2));
        // expand bounding box by threshold amount to ensure finding intersects at segment end points
        treeIndex.query(bgi::intersects(bb), queryResults);

        for (std::size_t i1 : queryResults) {
            std::size_t j1 = nextWrappingIndex(i1, pline1);
            Vertex<T> const &p1v1 = pline1[i1];
            Vertex<T> const &p1v2 = pline1[j1];

            auto intrAtStartPt = [&](Point<T> const &intr) {
                return point::is_equal(vertex::get_pos(p1v1), intr) || point::is_equal(vertex::get_pos(p2v1), intr);
            };

            auto intrResult = intrPlineSegs(p1v1, p1v2, p2v1, p2v2);
            switch (intrResult.intrType) {
            case PlineSegIntrType::NoIntersect:
                break;
            case PlineSegIntrType::TangentIntersect:
            case PlineSegIntrType::OneIntersect:
                if (!intrAtStartPt(intrResult.point1)) {
                    intrs.emplace_back(i1, i2, intrResult.point1);
                }
                break;
            case PlineSegIntrType::TwoIntersects:
                if (!intrAtStartPt(intrResult.point1)) {
                    intrs.emplace_back(i1, i2, intrResult.point1);
                }
                if (!intrAtStartPt(intrResult.point2)) {
                    intrs.emplace_back(i1, i2, intrResult.point2);
                }
                break;
            case PlineSegIntrType::SegmentOverlap:
            case PlineSegIntrType::ArcOverlap:
                coincidentIntrs.emplace_back(i1, i2, intrResult.point1, intrResult.point2);
                if (point::is_equal(vertex::get_pos(p1v1), intrResult.point1) ||
                    point::is_equal(vertex::get_pos(p1v1), intrResult.point2)) {
                    possibleDuplicates.insert({prevWrappingIndex(i1, pline1), i2});
                }
                if (point::is_equal(vertex::get_pos(p2v1), intrResult.point1) ||
                    point::is_equal(vertex::get_pos(p2v1), intrResult.point2)) {
                    possibleDuplicates.insert({i1, prevWrappingIndex(i2, pline2)});
                }
                break;
            }
        }

        // visit all indexes
        return true;
    };

    pline2.visitSegIndices(pline2SegVisitor);

    // remove duplicate points caused by the coincident intersect definition
    intrs.erase(std::remove_if(intrs.begin(), intrs.end(),
                               [&](auto const &intr) {
                                   auto found = possibleDuplicates.find({intr.sIndex1, intr.sIndex2});
                                   if (found == possibleDuplicates.end()) {
                                       return false;
                                   }

                                   auto const &endPt1 =
                                       vertex::get_pos(pline1[nextWrappingIndex(intr.sIndex1, pline1)]);
                                   if (fuzzyEqual(intr.pos, endPt1)) {
                                       return true;
                                   }

                                   auto const &endPt2 =
                                       vertex::get_pos(pline2[nextWrappingIndex(intr.sIndex2, pline2)]);
                                   return point::is_equal(vertex::get_pos(intr.pos), endPt2);
                               }),
                intrs.end());
}

template <typename T>
void offsetCircleIntersectsWithPline(Polyline<T> const &pline, jereal offset,
                                     Point<T> const &circleCenter,
                                     bmodel::BBoxRTree<T> const &treeIndex,
                                     std::vector<std::pair<std::size_t, Point<T>>> &output) {

    const T circleRadius = jemath::jeabs(offset);

    std::vector<std::size_t> queryResults;

    treeIndex.query(px(circleCenter) - circleRadius, py(circleCenter) - circleRadius,
                       px(circleCenter) + circleRadius, py(circleCenter) + circleRadius,
                       queryResults);

    auto validLineSegIntersect = [](T t) {
        return !falseIntersect(t) && jemath::jeabs(t) > J_V0;
    };

    auto validArcSegIntersect = [](Point<T> const &arcCenter, Point<T> const &arcStart,
                                   Point<T> const &arcEnd, T bulge,
                                   Point<T> const &intrPoint) {
        return (point::is_equal(arcStart, intrPoint)) && pointWithinArcSweepAngle(arcCenter, arcStart, arcEnd, bulge, intrPoint);
    };

    for (std::size_t sIndex : queryResults) {
        Vertex<T> const &v1 = pline[sIndex];
        Vertex<T> const &v2 = pline[sIndex + 1];
        if (jemath::fuzzy_iszero(vertex::get_angle(v1))) {
            IntrLineSeg2Circle2Result intrResult =
                intrLineSeg2Circle2(vertex::get_pos(v1), vertex::get_pos(v2), circleRadius, circleCenter);
            if (intrResult.numIntersects == 0) {
                continue;
            } else if (intrResult.numIntersects == 1) {
                if (validLineSegIntersect(intrResult.t0)) {
                    output.emplace_back(sIndex, pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t0));
                }
            } else {
                if (validLineSegIntersect(intrResult.t0)) {
                    output.emplace_back(sIndex, pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t0));
                }
                if (validLineSegIntersect(intrResult.t1)) {
                    output.emplace_back(sIndex, pointFromParametric(vertex::get_pos(v1), vertex::get_pos(v2), intrResult.t1));
                }
            }
        } else {
            auto arc = arcRadiusAndCenter(v1, v2);
            IntrCircle2Circle2Result intrResult =
                intrCircle2Circle2(arc.radius, arc.center, circleRadius, circleCenter);
            switch (intrResult.intrType) {
            case Circle2Circle2IntrType::NoIntersect:
                break;
            case Circle2Circle2IntrType::OneIntersect:
                if (validArcSegIntersect(arc.center, vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_angle(v1), intrResult.point1)) {
                    output.emplace_back(sIndex, intrResult.point1);
                }
                break;
            case Circle2Circle2IntrType::TwoIntersects:
                if (validArcSegIntersect(arc.center, vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_angle(v1), intrResult.point1)) {
                    output.emplace_back(sIndex, intrResult.point1);
                }
                if (validArcSegIntersect(arc.center, vertex::get_pos(v1), vertex::get_pos(v2), vertex::get_angle(v1), intrResult.point2)) {
                    output.emplace_back(sIndex, intrResult.point2);
                }
                break;
            case Circle2Circle2IntrType::Coincident:
                break;
            }
        }
    }
}

/// Function to test if a point is a valid distance from the original polyline.
template <typename T, std::size_t N>
bool pointValidForOffset(Polyline<T> const &pline, jereal offset,
                         bmodel::BBoxRTree<T> const &treeIndex,
                         Point<T> const &point, std::vector<std::size_t> &queryResults
                         ) {
    const jereal absOffset = jemath::jeabs(offset);
    const jereal minDist = absOffset * absOffset;

    bool pointValid = true;

    auto visitor = [&](std::size_t i) {
        std::size_t j = nextWrappingIndex(i, polyline::vertices(pline)); // index += 1
        auto closestPoint = closestPointOnSeg(pline[i], pline[j], point);
        jereal dist = distSquared(closestPoint, point);
        pointValid = dist > minDist;
        return pointValid;
    };

    treeIndex.query(bgi::intersects(get_bbox(px(point) - absOffset, py(point) - absOffset,
                                             px(point) + absOffset, py(point) + absOffset))
                    , queryResults);
    return pointValid;
}

template<typename T>
std::vector<PoylineOffsetSegment<T>> createUntrimmedOffsetSegments(Polyline<T> const &polyline,
                                                                jereal offset) {
    std::size_t segmentCount = polyline::verteics_count(polyline) - 1;

    std::vector<PoylineOffsetSegment<T>> result;
    result.reserve(segmentCount);

    auto lineVisitor = [&](Vertex<T> const &v1, Vertex<T> const &v2) {
        result.emplace_back();
        PoylineOffsetSegment<T> &seg = result.back();
        seg.collapsedArc = false;
        seg.origV2Pos = vertex::get_pos(v2);
        Point<T> edge = point_sub(vertex::get_pos(v2), vertex::get_pos(v1));
        Point<T> offsetV = unitPerp(edge) * offset;
        set_pos(seg.v1, point_add(vertex::get_pos(v1), offsetV));
        set_pos(seg.v2, point_add(vertex::get_pos(v2), offsetV));
    };

    auto arcVisitor = [&](Vertex<T> const &v1, Vertex<T> const &v2) {
        auto arc = arcRadiusAndCenter(v1, v2);
        T offs = vertex::get_angle(v2) < J_V0 ? offset : -offset;
        T radiusAfterOffset = arc.radius + offs;
        Point<T> v1ToCenter = point_sub(vertex::get_pos(v1), arc.center);
        normalize(v1ToCenter);
        Point<T> v2ToCenter = point_sub(vertex::get_pos(v2), arc.center);
        normalize(v2ToCenter);

        result.emplace_back();
        PoylineOffsetSegment<T> &seg = result.back();
        seg.origV2Pos = vertex::get_pos(v2);
        set_pos(seg.v1, point_add(point_mul(v1ToCenter, offs), vertex::get_pos(v1)));
        set_pos(seg.v2, point_add(point_mul(v2ToCenter, offs), vertex::get_pos(v2)));
        set_angle(seg.v2, vertex::get_angle(v2));

        if (jemath::fuzzy_iszero(vertex::get_angle(v1))) {
            // collapsed arc, offset arc start and end points towards arc center and turn into line
            // handles case where offset vertexes are equal and simplifies path for clipping algorithm
            seg.collapsedArc = true;
        } else {
            seg.collapsedArc = false;
            set_angle(seg.v2, vertex::get_angle(v2));
        }
    };

    auto offsetVisitor = [&](Vertex<T> const &v1, Vertex<T> const &v2) {
        if (jemath::fuzzy_iszero(vertex::get_angle(v2))) {
            lineVisitor(v1, v2);
        } else {
            arcVisitor(v1, v2);
        }
    };

    for (std::size_t i = 1; i < polyline::verteics_count(polyline); ++i) {
        offsetVisitor(polyline[i - 1], polyline[i]);
    }

    return result;
}

template<typename T>
Polyline<T> createRawOffsetPolyline(Polyline<T> const &polyline, jereal offset) {

    Polyline<T> result;
    if (polyline::verteics_count(result) < 2) {
        return result;
    }

    std::vector<PoylineOffsetSegment<T>> rawOffsets = createUntrimmedOffsetSegments(polyline, offset);
    if (rawOffsets.size() == 0) {
        return result;
    }

    // detect single collapsed arc segment (this may be removed in the future if invalid segments are
    // tracked in join functions to be pruned at slice creation)
    if (rawOffsets.size() == 1 && rawOffsets[0].collapsedArc) {
        return result;
    }

    polyline::vertices(result).reserve(polyline::verteics_count(polyline));

    const bool connectionArcsAreCW = offset > static_cast<T>(0);

    auto joinResultVisitor = [connectionArcsAreCW](PoylineOffsetSegment<T> const &s1,
                                                   PoylineOffsetSegment<T> const &s2,
                                                   Polyline<T> &result) {
        const bool s1IsLine = (jemath::fuzzy_iszero(vertex::get_angle(s1.v1)));
        const bool s2IsLine = (jemath::fuzzy_iszero(vertex::get_angle(s2.v1)));
        if (s1IsLine && s2IsLine) {
            lineToLineJoin(s1, s2, connectionArcsAreCW, result);
        } else if (s1IsLine) {
            lineToArcJoin(s1, s2, connectionArcsAreCW, result);
        } else if (s2IsLine) {
            arcToLineJoin(s1, s2, connectionArcsAreCW, result);
        } else {
            arcToArcJoin(s1, s2, connectionArcsAreCW, result);
        }
    };

    polyline::pushback_vertex(result, rawOffsets[0].v1);

    // join first two segments and determine if first vertex was replaced (to know how to handle last
    // two segment joins for closed polyline)
    if (rawOffsets.size() > 1) {
        auto const &seg01 = rawOffsets[0];
        auto const &seg12 = rawOffsets[1];
        joinResultVisitor(seg01, seg12, result);
    }
    const bool firstVertexReplaced = (polyline::verteics_count(result) == 1);

    for (std::size_t i = 2; i < rawOffsets.size(); ++i) {
        const auto &seg1 = rawOffsets[i - 1];
        const auto &seg2 = rawOffsets[i];
        joinResultVisitor(seg1, seg2, result);
    }

    if (polyline::is_closed(polyline) && polyline::verteics_count(result) > 1) {
        // joining segments at vertex indexes (n, 0) and (0, 1)
        const auto &s1 = rawOffsets.back();
        const auto &s2 = rawOffsets[0];

        // temp polyline to capture results of joining (to avoid mutating result)
        Polyline<T> closingPartResult;
        polyline::pushback_vertex(closingPartResult, polyline::last_vertex(result));
        joinResultVisitor(s1, s2, closingPartResult);

        // update last vertexes
        polyline::last_vertex(result) = closingPartResult[0];
        for (std::size_t i = 1; i < polyline::verteics_count(closingPartResult); ++i) {
            polyline::pushback_vertex(result, closingPartResult[i]);
        }
        polyline::popback_vertex(result);

        // update first vertex (only if it has not already been updated/replaced)
        if (!firstVertexReplaced) {
            const Point<T> &updatedFirstPos = vertex::get_pos(polyline::last_vertex(closingPartResult));
            if (jemath::fuzzy_iszero(vertex::get_angle(result[0]))) {
                // just update position
                set_pos((result[0]), updatedFirstPos);
            } else if (polyline::verteics_count(result) > 1) {
                // update position and bulge
                const auto arc = arcRadiusAndCenter(result[0], result[1]);
                const jereal a1 = calc_angle(arc.center, updatedFirstPos);
                const jereal a2 = calc_angle(arc.center, vertex::get_pos(result[1]));
                const jereal updatedTheta = deltaAngle(a1, a2);
                if ((updatedTheta < J_V0 && vertex::get_angle(result[0]) > J_V0) ||
                    (updatedTheta > J_V0 && vertex::get_angle(result[0]) < J_V0)) {
                    // first vertex not valid, just update its position to be removed later
                    set_pos(result[0], updatedFirstPos);
                } else {
                    // update position and bulge
                    set_pos(result[0], updatedFirstPos);
                    set_angle(result[0], updatedTheta);
                }
            }
        }

        // must do final singularity prune between first and second vertex after joining curves (n, 0)
        // and (0, 1)
        if (polyline::verteics_count(result) > 1) {
            if (set_pos(result[0], vertex::get_pos(result[1]))) {
                polyline::remove_vertex(result, 0);
            }
        }
    } else {
        addOrReplaceIfSamePos(result, rawOffsets.back().v2);
    }

    // if due to joining of segments we are left with only 1 vertex then return no raw offset (empty
    // polyline)
    if (polyline::verteics_count(result) == 1) {
        result.clear();
    }

    return result;
}

/// Represents an open polyline slice of the raw offset polyline.
template<typename T>
struct OpenPolylineSlice {
    std::size_t intrStartIndex;
    Polyline<T> polyline;
    OpenPolylineSlice() = default;
    OpenPolylineSlice(std::size_t sIndex, Polyline<T> slice)
        : intrStartIndex(sIndex), polyline(std::move(slice)) {}
};

template <typename T>
std::vector<OpenPolylineSlice<T>> slicesFromRawOffset(Polyline<T> const &originalPline,
                                                   Polyline<T> const &rawOffsetPline,
                                                   jereal offset) {

    std::vector<OpenPolylineSlice<T>> result;
    if (polyline::verteics_count(rawOffsetPline) < 2) {
        return result;
    }

    bmodel::BBoxRTree<T> origPlineSpatialIndex = rtree::create_box_rtree(polyline::get_edges(originalPline));
    bmodel::BBoxRTree<T> rawOffsetPlineSpatialIndex = rtree::create_box_rtree(polyline::get_edges(rawOffsetPline));

    std::vector<PlineIntersect<T>> selfIntersects;
    allSelfIntersects(rawOffsetPline, selfIntersects, rawOffsetPlineSpatialIndex);

    std::vector<std::size_t> queryResults;
    if (selfIntersects.size() == 0) {
        if (!pointValidForOffset(originalPline, offset, origPlineSpatialIndex, vertex::get_pos(rawOffsetPline[0]),
                                 queryResults)) {
            return result;
        }
        // copy and convert raw offset into open polyline
        result.emplace_back(std::numeric_limits<std::size_t>::max(), rawOffsetPline);
        polyline::pushback_vertex(result.back().polyline, rawOffsetPline[0]);
        set_angle(polyline::last_vertex(result.back().polyline), J_V0);
        return result;
    }

    // using unordered_map rather than map for performance (as is used in
    // dualSliceAtIntersectsForOffset) since all slices will stitch together to form closed
    // loops/polylines so later when slices are stitched together the order that slices are visited
    // does not matter
    std::unordered_map<std::size_t, std::vector<Point<T>>> intersectsLookup;
    intersectsLookup.reserve(2 * selfIntersects.size());

    for (PlineIntersect<T> const &si : selfIntersects) {
        intersectsLookup[si.sIndex1].push_back(si.pos);
        intersectsLookup[si.sIndex2].push_back(si.pos);
    }

    // sort intersects by distance from start vertex
    for (auto &kvp : intersectsLookup) {
        Point<T> startPos = vertex::get_pos(rawOffsetPline[kvp.first]);
        auto cmp = [&](Point<T> const &si1, Point<T> const &si2) {
            return distSquared(si1, startPos) < distSquared(si2, startPos);
        };
        std::sort(kvp.second.begin(), kvp.second.end(), cmp);
    }

    auto intersectsOrigPline = [&](Vertex<T> const &v1, Vertex<T> const &v2) {
        bmodel::BBox<T> approxBB = rect::to_bbox(get_bound_rect(v1, v2));
        bool hasIntersect = false;
        auto visitor = [&](std::size_t i) {
            std::size_t j = nextWrappingIndex(i, originalPline);
            IntrPlineSegsResult<T> intrResult =
                intrPlineSegs(v1, v2, originalPline[i], originalPline[j]);
            hasIntersect = intrResult.intrType != PlineSegIntrType::NoIntersect;
            return !hasIntersect;
        };

        origPlineSpatialIndex.query(bgi::intersects(approxBB), queryResults);

        return hasIntersect;
    };

    for (auto const &kvp : intersectsLookup) {
        // start index for the slice we're about to build
        std::size_t sIndex = kvp.first;
        // self intersect list for this start index
        std::vector<Point<T>> const &siList = kvp.second;

        const auto &startVertex = rawOffsetPline[sIndex];
        std::size_t nextIndex = nextWrappingIndex(sIndex, rawOffsetPline);
        const auto &endVertex = rawOffsetPline[nextIndex];

        if (siList.size() != 1) {
            // build all the segments between the N intersects in siList (N > 1), skipping the first
            // segment (to be processed at the end)
            SplitResult<T> firstSplit = splitAtPoint(startVertex, endVertex, siList[0]);
            auto prevVertex = firstSplit.splitVertex;
            for (std::size_t i = 1; i < siList.size(); ++i) {
                SplitResult<T> split = splitAtPoint(prevVertex, endVertex, siList[i]);
                // update prevVertex for next loop iteration
                prevVertex = split.splitVertex;
                // skip if they're ontop of each other
                if (point::is_equal(vertex::get_pos(split.updatedStart), vertex::get_pos(split.splitVertex))) {
                    continue;
                }

                // test start point
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   vertex::get_pos(split.updatedStart), queryResults)) {
                    continue;
                }

                // test end point
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   vertex::get_pos(split.splitVertex), queryResults)) {
                    continue;
                }

                // test mid point
                auto midpoint = segMidpoint(split.updatedStart, split.splitVertex);
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, midpoint,
                                                   queryResults)) {
                    continue;
                }

                // test intersection with original polyline
                if (intersectsOrigPline(split.updatedStart, split.splitVertex)) {
                    continue;
                }

                result.emplace_back();
                result.back().intrStartIndex = sIndex;
                polyline::pushback_vertex(result.back().pline, split.updatedStart);
                polyline::pushback_vertex(result.back().pline, split.splitVertex);
            }
        }

        // build the segment between the last intersect in siList and the next intersect found

        // check that the first point is valid
        if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, siList.back(),
                                           queryResults)) {
            continue;
        }

        SplitResult<T> split = splitAtPoint(startVertex, endVertex, siList.back());
        Polyline<T> currSlice;
        polyline::pushback_vertex(currSlice, split.splitVertex);

        std::size_t index = nextIndex;
        bool isValidPline = true;
        std::size_t loopCount = 0;
        const std::size_t maxLoopCount = rawOffsetPline.size();
        while (true) {
            if (loopCount++ > maxLoopCount) {
                // break to avoid infinite loop
                break;
            }
            // check that vertex point is valid
            if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                               vertex::get_pos(rawOffsetPline[index]), queryResults)) {
                isValidPline = false;
                break;
            }

            // check that the segment does not intersect original polyline
            if (intersectsOrigPline(polyline::last_vertex(currSlice), rawOffsetPline[index])) {
                isValidPline = false;
                break;
            }

            // add vertex
            addOrReplaceIfSamePos(currSlice, rawOffsetPline[index]);

            // check if segment that starts at vertex we just added has an intersect
            auto nextIntr = intersectsLookup.find(index);
            if (nextIntr != intersectsLookup.end()) {
                // there is an intersect, slice is done, check if final segment is valid

                // check intersect pos is valid (which will also be end vertex position)
                Point<T> const &intersectPos = nextIntr->second[0];
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   intersectPos, queryResults)) {
                    isValidPline = false;
                    break;
                }

                std::size_t nextIndex = nextWrappingIndex(index, rawOffsetPline);
                SplitResult<T> split =
                    splitAtPoint(polyline::last_vertex(currSlice), rawOffsetPline[nextIndex], intersectPos);

                Vertex<T> sliceEndVertex = Vertex<T>(intersectPos, T(0));
                // check mid point is valid
                Point<T> mp = segMidpoint(split.updatedStart, sliceEndVertex);
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, mp,
                                                   queryResults)) {
                    isValidPline = false;
                    break;
                }

                // trim last added vertex and add final intersect position
                polyline::last_vertex(currSlice) = split.updatedStart;
                addOrReplaceIfSamePos(currSlice, sliceEndVertex);

                break;
            }
            // else there is not an intersect, increment index and continue
            index = nextWrappingIndex(index, rawOffsetPline);
        }

        isValidPline = isValidPline && currSlice.size() > 1;

        if (isValidPline && point::is_equal(vertex::get_pos(currSlice[0]), vertex::get_pos(polyline::last_vertex(currSlice)))) {
            // discard very short slice loops (invalid loops may arise due to valid offset distance
            // thresholding)
            isValidPline = getPathLength(currSlice) > T(1e-2);
        }

        if (isValidPline) {
            result.emplace_back(sIndex, std::move(currSlice));
        }
    }

    return result;
}

/// Slices a raw offset polyline at all of its self intersects and intersects with its dual.
template <typename T>
std::vector<OpenPolylineSlice<T>>
dualSliceAtIntersectsForOffset(Polyline<T> const &originalPline,
                               Polyline<T> const &rawOffsetPline,
                               Polyline<T> const &dualRawOffsetPline, jereal offset) {
    std::vector<OpenPolylineSlice<T>> result;
    if (polyline::verteics_count(rawOffsetPline) < 2) {
        return result;
    }

    bmodel::BBoxRTree<T> origPlineSpatialIndex = rtree::create_box_rtree(polyline::get_edges(originalPline));
    bmodel::BBoxRTree<T> rawOffsetPlineSpatialIndex = rtree::create_box_rtree(polyline::get_edges(rawOffsetPline));

    std::vector<PlineIntersect<T>> selfIntersects;
    allSelfIntersects(rawOffsetPline, selfIntersects, rawOffsetPlineSpatialIndex);

    PlineIntersectsResult<T> dualIntersects;
    findIntersects(rawOffsetPline, dualRawOffsetPline, rawOffsetPlineSpatialIndex, dualIntersects);

    // using map rather than unordered map since we want to construct the slices in vertex index order
    // and we do so by looping through all intersects (required later when slices are stitched
    // together, because slices may not all form closed loops/polylines so must go in order of
    // indexes to ensure longest sitched results are formed)
    std::map<std::size_t, std::vector<Point<T>>> intersectsLookup;

    if (!polyline::is_closed(originalPline)) {
        // find intersects between circles generated at original open polyline end points and raw offset
        // polyline
        std::vector<std::pair<std::size_t, Point<T>>> intersects;
        internal::offsetCircleIntersectsWithPline(rawOffsetPline, offset, vertex::get_pos(originalPline[0]),
                                                  rawOffsetPlineSpatialIndex, intersects);
        internal::offsetCircleIntersectsWithPline(rawOffsetPline, offset,
                                                  vertex::get_pos(polyline::last_vertex(originalPline)),
                                                  rawOffsetPlineSpatialIndex, intersects);
        for (auto const &pair : intersects) {
            intersectsLookup[pair.first].push_back(pair.second);
        }
    }

    for (PlineIntersect<T> const &si : selfIntersects) {
        intersectsLookup[si.sIndex1].push_back(si.pos);
        intersectsLookup[si.sIndex2].push_back(si.pos);
    }

    for (PlineIntersect<T> const &intr : dualIntersects.intersects) {
        intersectsLookup[intr.sIndex1].push_back(intr.pos);
    }

    for (PlineCoincidentIntersect<T> const &intr : dualIntersects.coincidentIntersects) {
        intersectsLookup[intr.sIndex1].push_back(intr.point1);
        intersectsLookup[intr.sIndex1].push_back(intr.point2);
    }

    std::vector<std::size_t> queryResults;
    queryResults.reserve(8);
    if (intersectsLookup.size() == 0) {
        if (!pointValidForOffset(originalPline, offset, origPlineSpatialIndex, vertex::get_pos(rawOffsetPline[0]),
                                 queryResults)) {
            return result;
        }
        // copy and convert raw offset into open polyline
        result.emplace_back(std::numeric_limits<std::size_t>::max(), rawOffsetPline);
        if (polyline::is_closed(originalPline)) {
            polyline::pushback_vertex(result.back().polyline, rawOffsetPline[0]);
            set_angle(polyline::last_vertex(result.back().polyline), J_V0);
        }
        return result;
    }

    // sort intersects by distance from start vertex
    for (auto &kvp : intersectsLookup) {
        Point<T> startPos = vertex::get_pos(rawOffsetPline[kvp.first]);
        auto cmp = [&](Point<T> const &si1, Point<T> const &si2) {
            return distSquared(si1, startPos) < distSquared(si2, startPos);
        };
        std::sort(kvp.second.begin(), kvp.second.end(), cmp);
    }

    auto intersectsOrigPline = [&](Vertex<T> const &v1, Vertex<T> const &v2) {
        bmodel::BBox<T> approxBB = rect::to_bbox(get_bound_rect(v1, v2));
        bool intersects = false;
        auto visitor = [&](std::size_t i) {
            std::size_t j = nextWrappingIndex(i, originalPline);
            IntrPlineSegsResult<T> intrResult =
                intrPlineSegs(v1, v2, originalPline[i], originalPline[j]);
            intersects = intrResult.intrType != PlineSegIntrType::NoIntersect;
            return !intersects;
        };

        origPlineSpatialIndex.query(bgi::intersects(approxBB), queryResults);

        return intersects;
    };

    if (!polyline::is_closed(originalPline)) {
        // build first open polyline that ends at the first intersect since we will not wrap back to
        // capture it as in the case of a closed polyline
        Polyline<T> firstSlice;
        std::size_t index = 0;
        std::size_t loopCount = 0;
        const std::size_t maxLoopCount = polyline::verteics_count(rawOffsetPline);
        while (true) {
            if (loopCount++ > maxLoopCount) {
                // break to avoid infinite loop
                break;
            }
            auto iter = intersectsLookup.find(index);
            if (iter == intersectsLookup.end()) {
                // no intersect found, test segment will be valid before adding the vertex
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   vertex::get_pos(rawOffsetPline[index]), queryResults)) {
                    break;
                }

                // index check (only test segment if we're not adding the first vertex)
                if (index != 0 && intersectsOrigPline(polyline::last_vertex(firstSlice), rawOffsetPline[index])) {
                    break;
                }

                addOrReplaceIfSamePos(firstSlice, rawOffsetPline[index]);
            } else {
                // intersect found, test segment will be valid before finishing first open polyline
                Point<T> const &intersectPos = iter->second[0];
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   intersectPos, queryResults)) {
                    break;
                }

                SplitResult<T> split =
                    splitAtPoint(rawOffsetPline[index], rawOffsetPline[index + 1], intersectPos);

                Vertex<T> sliceEndVertex = Vertex<T>(intersectPos, T(0));
                auto midpoint = segMidpoint(split.updatedStart, sliceEndVertex);
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, midpoint,
                                                   queryResults)) {
                    break;
                }

                if (intersectsOrigPline(split.updatedStart, sliceEndVertex)) {
                    break;
                }

                addOrReplaceIfSamePos(firstSlice, split.updatedStart);
                addOrReplaceIfSamePos(firstSlice, sliceEndVertex);
                result.emplace_back(0, std::move(firstSlice));
                break;
            }

            index += 1;
        }
    }

    for (auto const &kvp : intersectsLookup) {
        // start index for the slice we're about to build
        std::size_t sIndex = kvp.first;
        // self intersect list for this start index
        std::vector<Point<T>> const &siList = kvp.second;

        const auto &startVertex = rawOffsetPline[sIndex];
        std::size_t nextIndex = nextWrappingIndex(sIndex, rawOffsetPline);
        const auto &endVertex = rawOffsetPline[nextIndex];

        if (siList.size() != 1) {
            // build all the segments between the N intersects in siList (N > 1), skipping the first
            // segment (to be processed at the end)
            SplitResult<T> firstSplit = splitAtPoint(startVertex, endVertex, siList[0]);
            auto prevVertex = firstSplit.splitVertex;
            for (std::size_t i = 1; i < siList.size(); ++i) {
                SplitResult<T> split = splitAtPoint(prevVertex, endVertex, siList[i]);
                // update prevVertex for next loop iteration
                prevVertex = split.splitVertex;
                // skip if they're ontop of each other
                if (point::is_equal(vertex::get_pos(split.updatedStart), vertex::get_pos(split.splitVertex),
                               realPrecision<T>())) {
                    continue;
                }

                // test start point
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   vertex::get_pos(split.updatedStart), queryResults)) {
                    continue;
                }

                // test end point
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   vertex::get_pos(split.splitVertex), queryResults)) {
                    continue;
                }

                // test mid point
                auto midpoint = segMidpoint(split.updatedStart, split.splitVertex);
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, midpoint,
                                                   queryResults)) {
                    continue;
                }

                // test intersection with original polyline
                if (intersectsOrigPline(split.updatedStart, split.splitVertex)) {
                    continue;
                }

                result.emplace_back();
                result.back().intrStartIndex = sIndex;
                polyline::pushback_vertex(result.back().polyline, split.updatedStart);
                polyline::pushback_vertex(result.back().polyline, split.splitVertex);
            }
        }

        // build the segment between the last intersect in siList and the next intersect found

        // check that the first point is valid
        if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, siList.back(),
                                           queryResults)) {
            continue;
        }

        SplitResult<T> split = splitAtPoint(startVertex, endVertex, siList.back());
        Polyline<T> currSlice;
        polyline::pushback_vertex(currSlice, split.splitVertex);

        std::size_t index = nextIndex;
        bool isValidPline = true;
        std::size_t loopCount = 0;
        const std::size_t maxLoopCount = rawOffsetPline.size();
        while (true) {
            if (loopCount++ > maxLoopCount) {
                // break to avoid infinite loop
                break;
            }
            // check that vertex point is valid
            if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                               vertex::get_pos(rawOffsetPline[index]), queryResults)) {
                isValidPline = false;
                break;
            }

            // check that the segment does not intersect original polyline
            if (intersectsOrigPline(polyline::last_vertex(currSlice), rawOffsetPline[index])) {
                isValidPline = false;
                break;
            }

            // add vertex
            addOrReplaceIfSamePos(currSlice, rawOffsetPline[index]);

            // check if segment that starts at vertex we just added has an intersect
            auto nextIntr = intersectsLookup.find(index);
            if (nextIntr != intersectsLookup.end()) {
                // there is an intersect, slice is done, check if final segment is valid

                // check intersect pos is valid (which will also be end vertex position)
                Point<T> const &intersectPos = nextIntr->second[0];
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex,
                                                   intersectPos, queryResults)) {
                    isValidPline = false;
                    break;
                }

                std::size_t nextIndex = nextWrappingIndex(index, rawOffsetPline);
                SplitResult<T> split =
                    splitAtPoint(polyline::last_vertex(currSlice), rawOffsetPline[nextIndex], intersectPos);

                Vertex<T> sliceEndVertex = Vertex<T>(intersectPos, T(0));
                // check mid point is valid
                Point<T> mp = segMidpoint(split.updatedStart, sliceEndVertex);
                if (!internal::pointValidForOffset(originalPline, offset, origPlineSpatialIndex, mp,
                                                   queryResults)) {
                    isValidPline = false;
                    break;
                }

                // trim last added vertex and add final intersect position
                polyline::last_vertex(currSlice) = split.updatedStart;
                addOrReplaceIfSamePos(currSlice, sliceEndVertex);

                break;
            }
            // else there is not an intersect, increment index and continue
            if (index == polyline::verteics_count(rawOffsetPline) - 1) {
                if (polyline::is_closed(originalPline)) {
                    // wrap index
                    index = 0;
                } else {
                    // open polyline, we're done
                    break;
                }
            } else {
                index += 1;
            }
        }

        if (isValidPline && polyline::verteics_count(currSlice) > 1) {
            result.emplace_back(sIndex, std::move(currSlice));
        }
    }

    return result;
}

/// Stitches raw offset polyline slices together, discarding any that are not valid.
template <typename T>
std::vector<Polyline<T>>
stitchOffsetSlicesTogether(std::vector<OpenPolylineSlice<T>> const &slices, bool closedPolyline,
                           std::size_t origMaxIndex,
                           T joinThreshold = jemath::Epsilon<T>()) {
    std::vector<Polyline<T>> result;
    if (slices.size() == 0) {
        return result;
    }

    if (slices.size() == 1) {
        result.emplace_back(slices[0].polyline);
        if (closedPolyline &&
            point::is_equal(vertex::get_pos(result[0][0]), vertex::get_pos(polyline::last_vertex(result[0])), joinThreshold)) {
            polyline::popback_vertex(result[0]);
        }

        return result;
    }

    // load spatial index with all start points
    bmodel::BBoxRTree<T> treeIndex;

    size_t n = 0;
    for (const auto &slice : slices) {
        auto const &point = vertex::get_pos(slice.polyline[0]);
        treeIndex.insert(std::make_pair(get_bbox(px(point) - joinThreshold, py(point) - joinThreshold,
                                                 px(point) + joinThreshold, py(point) + joinThreshold), n));
        n++;
    }

    std::vector<bool> visitedIndexes(slices.size(), false);
    std::vector<std::size_t> queryResults;
    for (std::size_t i = 0; i < slices.size(); ++i) {
        if (visitedIndexes[i]) {
            continue;
        }

        visitedIndexes[i] = true;

        Polyline<T> currPline;
        std::size_t currIndex = i;
        auto const &initialStartPoint = vertex::get_pos(slices[i].polyline[0]);
        std::size_t loopCount = 0;
        const std::size_t maxLoopCount = slices.size();
        while (true) {
            if (loopCount++ > maxLoopCount) {
                // break to avoid infinite loop
                break;
            }
            const std::size_t currLoopStartIndex = slices[currIndex].intrStartIndex;
            auto const &currSlice = slices[currIndex].polyline;
            auto const &currEndPoint = vertex::get_pos(polyline::last_vertex(currSlice));
            polyline::vertices(currPline).insert(polyline::vertices(currPline).end(), polyline::vertices(currSlice).begin(),
                                       polyline::vertices(currSlice).end());
            queryResults.clear();
            Point<T> leftbottom{};
            Point<T> righttop{};
            treeIndex.query(bgi::intersects(get_bbox(px(currEndPoint) - joinThreshold, py(currEndPoint) - joinThreshold,
                                                     px(currEndPoint) + joinThreshold, py(currEndPoint) + joinThreshold))
                            , queryResults);

            queryResults.erase(std::remove_if(queryResults.begin(), queryResults.end(),
                                              [&](std::size_t index) { return visitedIndexes[index]; }),
                               queryResults.end());

            auto indexDistAndEqualInitial = [&](std::size_t index) {
                auto const &slice = slices[index];
                std::size_t indexDist;
                if (currLoopStartIndex <= slice.intrStartIndex) {
                    indexDist = slice.intrStartIndex - currLoopStartIndex;
                } else {
                    // forward wrapping distance (distance to end + distance to index)
                    indexDist = origMaxIndex - currLoopStartIndex + slice.intrStartIndex;
                }

                bool equalToInitial = (point::is_equal(vertex::get_pos(polyline::last_vertex(slice.pline)), initialStartPoint));

                return std::make_pair(indexDist, equalToInitial);
            };

            std::sort(queryResults.begin(), queryResults.end(),
                      [&](std::size_t index1, std::size_t index2) {
                          auto distAndEqualInitial1 = indexDistAndEqualInitial(index1);
                          auto distAndEqualInitial2 = indexDistAndEqualInitial(index2);
                          if (distAndEqualInitial1.first == distAndEqualInitial2.first) {
                              // index distances are equal, compare on position being equal to initial start
                              // (testing index1 < index2, we want the longest closed loop possible)
                              return distAndEqualInitial1.second < distAndEqualInitial2.second;
                          }

                          return distAndEqualInitial1.first < distAndEqualInitial2.first;
                      });

            if (queryResults.size() == 0) {
                // we're done
                if (polyline::verteics_count(currPline) > 1) {
                    if (closedPolyline && point::is_equal(
                            vertex::get_pos(currPline[0]), vertex::get_pos(polyline::last_vertex(currPline)))) {
                        polyline::popback_vertex(currPline);
                    }
                    result.emplace_back(std::move(currPline));
                }
                break;
            }

            // else continue stitching
            visitedIndexes[queryResults[0]] = true;
            polyline::popback_vertex(currPline);
            currIndex = queryResults[0];
        }
    }

    return result;
}
}


namespace offset {

template<typename T>
Polygon<T> polygon_offset(Polygon<T> const &poly, jereal offset) {
    using namespace internal;
    Polygon<T> result;
    for(auto path: poly)
    {
        if (polyline::verteics_count(path) < 2) {
            polygon::add_subpath(result, path);
            continue;
        }
        auto rawOffset = createRawOffsetPolyline(path, offset);
        if (polyline::is_closed(path)) {
            auto slices = slicesFromRawOffset(path, rawOffset, offset);
            polygon::add_subpaths(result, stitchOffsetSlicesTogether(
                                              slices, polyline::is_closed(path), polyline::verteics_count(rawOffset) - 1));
            continue;
        }

        // not closed polyline
        auto dualRawOffset = createRawOffsetPolyline(path, -offset);
        auto slices = dualSliceAtIntersectsForOffset(path, rawOffset, dualRawOffset, offset);
        polygon::add_subpaths(result, stitchOffsetSlicesTogether(
                                          slices, polyline::is_closed(path), polyline::verteics_count(rawOffset) - 1));
    }
    return result;
}

} // offset

} // jeo
#endif // OFFSET_HPP
